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Abstract 

We present systematic study of zero modes and gaps by introducing effects of anisotropy of 
hopping integrals for a tight-binding model on the honeycomb lattice in a magnetic field. The 
condition for the existence of zero modes is analytically derived. From the condition, it is found 
that a tiny anisotropy for graphene is sufficient to open a gap around zero energy in a magnetic 
field. This gap behaves as a non-perturbative and exponential form as a function of the magnetic 
field. The non-analytic behavior with respect to the magnetic field can be understood as tunneling 
effects between energy levels around two Dirac zero modes appearing in the honeycomb lattice, 
and an explicit form of the gap around zero energy is obtained by the WKB method near the 
merging point of these Dirac zero modes. Effects of the anisotropy for the honeycomb lattices 
with boundaries are also studied. The condition for the existence of zero energy edge states in a 
magnetic field is analytically derived. On the basis of the condition, it is recognized that anisotropy 
of the hopping integrals induces abrupt changes of the number of zero energy edge states, which 
depend on the shapes of the edges sensitively. 

PACS numbers: 71.70.Di, 73.43.-f, 81.05.Uw 
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FIG. 1: The honeycomb lattice. The hopping integrals of the horizontal bonds are t, and those for 
the other bonds are 1. A magnetic flux 27r<I> is applied through the unit hexagon. 



INTRODUCTION 



Recent experiments on grapheneQ, I, I, Q, Q have led to renewed interest in physical 
properties of electrons on the honeycomb lattice. Despite its simple structure, the hon- 
eycomb lattice provides non-trivial physical phenomena which can not be observed in the 
ordinary square lattice. Among them, much attention has been paid to its peculiar dis- 
persion. In the absence of a magnetic field, the honeycomb lattice has E = zero modes 
at the corners K and K' of the Brillouin zone. By treating these zero modes as Dirac 
fermions, the unconventional quantization of the Hall conductance observed for graphene 



was , although the full proper theoretical treatment of the Hall conduc- 

tance on the honeycomb lattice was made very recently^]. Moreover, when the system has 
a boundary, there are E = edge modes localized on the boundary. The existence of the 
E = edge modes depends on a choice of the boundary, and for zigzag and bearded edges 



there occur large density o 



flat band structures 
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states localized on these edges at the Fermi energy due to their 
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191 . 120(]. In addition, it is suggested 
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that the E = edge modes induce charge accumulation on these edges 

In this paper, we study properties of these E = zero modes in the presence of anisotropy 
of the hopping integrals in the honeycomb lattice. Recently, the anisotropy of the hopping 
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integrals was introduced by replacing one of the hopping integrals with a general value 



mm, 
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in order to investigate the unconventional quantum Hall effects 
on graphene (see Fig{T]). (For t 7^ 1, we have the anisotropic honeycomb lattice.) In Ref.j^, 
by using topological arguments, an algebraic expression of the quantum Hall conductance 
was obtained for almost all gaps including subband gaps, and it was shown that the uncon- 
ventional quantization of the Hall conductance in a weak magnetic field is realized for weak 
^ (0 < ^ ~ 1)) while only the conventional quantization is obtained for strong t {t > 2). 
Furthermore, for the graphene case {t = 1), the unconventional quantization was found to 
persist up to the Van Hove singularity [ol, Q. 

lopping parameters is also known to change the peculiar dispersion 



The anisotropy o: 



mentioned above 
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the 



25 



27l |. However, in the absence of a magnetic field, its influence 



is restrictive: Although a gap opens around zero energy for t > 2, there remain two E = 
zero modes in the Brillouin zone for < t < 2 [24]. Therefore, a large anisotropy is needed 
to change the zero mode structure. As well as the zero mode structure, that of zero energy 
edge states was shown to change by a large anisotropy [28 1. 

In this paper, it will be shown the situation is drastically changed in the presence of a 
magnetic field. We analytically derive the condition for the existence of zero modes in a 
magnetic flux 27r$ = 2Txp/q {p and q are mutually prime integers), and from the condition it 
is found that, in the limit of g — * cxd, zero modes exist only for < t < 1, but a gap around 
zero energy opens for t > 1. In other words, a small anisotropy t = l + e(0<e^l)is 
sufficient to open a gap in the presence of a weak magnetic field. 

For 1 < t < 2, the gap around zero energy in a weak magnetic field behaves as a non- 



perturbative and exponential form as a function of $. It will be shown that this 



Dehavior 



is naturally explained in terms of the spontaneous breaking of supersymmetry |29l . |30| . In 
particular, an explicit form of the gap around zero energy for t ~ 2 is obtained by the WKB 
method. At t = 2, the gap around zero energy in a weak magnetic field is found to make 
a transition from an exponential (non-perturbative) to a power-law (perturbative) behavior 
as a function of $, and for t > 2, the energy bands in a weak magnetic field show linear 
dependence on $. 

We will also show that the structure of E = edge states in the presence of a magnetic 
field is different from that in the absence of a magnetic field. The condition for the existence 
of zero energy edge states in a magnetic field is analytically derived, and it is found that 
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the anisotropy of the hopping integrals induces abrupt changes of the number of zero energy 
edge states, which also sensitively depend on shapes of the edges. 

The organization of this paper is as follows. In SecHIl we present our model. The 
condition for the existence of zero modes in a magnetic field is analytically derived in Sec JIIIl 
both from the secular equation and from the normalizability condition of states with zero 
energy. On the basis of the condition for the existence of zero modes, the energy spectrum 
near zero energy in a weak magnetic field is systematically examined in Sec JIVI In SecJVl zero 
energy edge states are analyzed, where crucial roles of the anisotropy of the hopping integrals 
are recognized again. Finally, we summarize our results and discuss possible experimental 
realization of anisotropy of the hopping integrals in Sec JVIl 

II. TIGHT-BINDING MODEL ON THE HONEYCOMB LATTICE IN A MAG- 
NETIC FIELD 

Let us consider the tight-binding model on the honeycomb lattice with nearest-neighbor 
hopping in a magnetic field as shown in Fig{Tl By denoting wave functions on two sublattices 
of the honeycomb lattice as ipn,m and 4>n,m, respectively, the tight-binding model is given by 

Elpn^m = 0n+l,m-l + e^*'^*"0ri+l,m+l + t4'n,m, 

E(f)n,m = 1pn-l,m+l + C"^^''*^""^ Vn-l,m-l + ttljn,m, (l) 

where a magnetic flux through the unit hexagon is given by 27r$. Here we have introduced 
anisotropy of the hopping integrals: The hopping integrals of the horizontal bonds are t, 
and those for the other bonds are 1. For simplicity, we neglect the spin degrees of freedom 
in the following. 

HI. THE CONDITION FOR THE EXISTENCE OF ZERO MODES 

For the isotropic case (t = 1), it was found that zero modes exist for all (rational) values 
of $^]. We now derive the condition for the existence of zero modes in the anisotropic 
case. 
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Before examining $ 7^ 0, let us first consider $ = 0[24^. For $ = 0, ([T]) gives 

E4'n,m = 0n+l,m-l + 0n+l,m+l + t(j)n,mi 
E(j)n,m = 4'n-l,m+l + 4'n-l,m-l 

From the Bloch's theorem, the wave functions are written as 41 1 



Substituting ([3]) into ([2]), we have 



V{k) =t + 2e'^^ cosfe 





V{k) 
V*{k) 
The eigenenergies E are given by 

E = ±\V{k)\ 



= ±y (t + 2 cos kx cos fcj^)^ + 4 sin^ /c^ cos^ /c^^. 
From ([6]) with ii^ = 0, we find two zero modes at 

fcj : (A:^, A;;+) = (^tt, cos-^ , fep- : (fc"", fej") = (^tt, - cos'^ , 
for < t < 2. By expanding and ky around and in ([7]), 

kx = k^^+Px, ky = kl^+py, (Ip^I, Ipj^I < 1), 

©(fc) is given by 

V±{p) = -itpx ± V4 - t2p^, 

where and T'-(p) are those near and respectively. From ([6]) and ([9]), 

dispersion relation of the Dirac zero mode is obtained: 



For t = 2, the two Dirac zero modes merge into a confluent point 

{k:,k*y) = in,0), 
and for t > 2, we have a gap around zero energy. 



A. Derivation from a secular equation 



Now we consider $ 7^ 0. We suppose that $ is a rational number, $ = p/q {p and q are 
mutually prime integers). Since Eq.([T]) has translational symmetry in the y-direction, the 
wave functions are written as 



^n,m = e^""^n, <Pn,m = e^""</'n, (12) 



and ([T]) becomes 



E0„ = (e*'^ + e-^'=-=^*-*("-i))^„_i + t^Pn. (13) 
By the gauge transformation ipn — >■ e'^^^Vn and 0„ — > e**^"0„, (fT3|) is rewritten as 

E(f)n = A:_i^„_i + tVn, (14) 

where An = 1 + exp[i(^i + 27r$ri)] with 6i = 2k. Since the spectrum is found to be 
invariant under the transformation 61^61 + ^ir/q, we can restrict the range of 9i to 
< 9i < 2,71 /q without loss of generality. Moreover, in (fT^ . {ipn,4>n) and {ipn+q,4>n+g) obey 
the same equation, thus from the Bloch's theorem, we have 

ipn+q = exp{iq92)ipn, 4>n+q = exp{iq92)4>n, (15) 

where ^2 satisfies < ^2 < 27r/g. Therefore f|T^ reduces to the eigenequation of a 2g x 2q 
matrix. In the secular equation of this, all non-constant terms containing less than q factors 
of e'^^ should cancel out each other since the eigenvalue has periodicity 27i/q with respect 
to 9i. From this property, it is found that the secular equation is written as the following 
form: 

FiE') + fi9,,92) = 0, (16) 
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where F{E'^) is a gth-order polynomial of E"^ with -F(O) = 0, and it is independent of (^i, 62)- 
The secular determinant for f|T4|) with E = determines /(^i,^2) as 

2 



( t Ai 

t A2 



det 



t A„_i 



^9 + (_l)'?+l2 COS (^^01 



t J 

Qn (1+ ^ 



ti + {-iy-^e"^'''l[An 



n=l 



(17) 



Here we have used 



n=l n=l n=l 



(18) 



which is derived from f{9i + 27r/g, 6*2) = f{Oi, 62). 
When t satisfies 

< t < 2^/'?, 



(19) 



the range of f {01,62) is < f {61,62) < {f + 2)^ and there exist two independent (^^i,6'2)'s 
with f{6i,62) = 0. From the secular equation (fT6l) . we have two E = modes at these 
(6'i,^^2)'s. On the other hand, if t satisfies 



t > 2^/^ 



(20) 



we have {f - 2)^ < f{6i, 62) < {f + If and there is no (^i, 62) with /(^i, 62) = 0. We have 
a gap around zero energy in this case. 

Here we note that the condition for the existence of zero modes for $ = 0, that is, 
< t < 2, is reproduced by f|T9|) with q = I. (When q = 1, reduces to that with $ = 0.) 



B. Derivation from the normalizability condition of states 



In Sec. IIII At we derived the condition for the existence of zero modes from the secular 
equation. Here, we re-derive it from the normalizability condition of states with zero energy. 
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Let us first consider (flU) with E = 



(21) 



Then for $ = p/q, ^ and ^ lead to 

q f 1 

t 



^Nq+l 



-in- 

\n=l 

v+i)7+«= (-7) (^n^nj^ 



1X9 



n I <P{N+l)q+l = 1 -7 I [1 + (-l)'^^^e*'^^']0{W+l)g+/. 



Nq+l 



1X9 

7 



(22) 



where / = 0,l,2,...,g — 1. Taking the absolute values of the both sides in (l22l) . we obtain 

,2 



'Nq+l I 



t9 



g^i + 1 

cos TT 

2 2 



\i^(N+l)q+l\ - ^ 



cosi— + ^vr 



l0(Af+l)g+«|5 



For t'? > 2, ([23]) gives 



\4'Nq+l\ < \(t){N+l)q+l\, \i'{N+l)q+l\ < IV'iVg+d- 



(23) 



(24) 



From (!24|) . we see that |0„| diverges for n ^ 00, and \ipn\ diverges for n — > — 00. Thus these 
states are not normalizable, and no relevant zero modes exist. We have a gap around E = Q 
in this case. On the other hand, for <2, (123!) gives 



^Nq+l\ 



\4>{N+1). 



q+l\} 



{N+l)q+l\ 



Nq+l I , 



(25) 



at 61 = ±|cos ^ f- + ^y'"^- Thus there exist two zero modes for t < 2^/'^. These results 
coincide with those of Sec. Illl Al 



IV. SPECTRUM NEAR ZERO ENERGY IN A WEAK MAGNETIC FIELD 



In this section, we examine the spectrum near zero energy in a weak magnetic field. 



Although some numerical study was presented in Ref.[23|, we perform detailed analytical 
study here. On the basis of the condition for the existence of zero modes obtained in the 
previous section, we consider the following four cases separately: 

A. < t < 1, where the condition (fT9|) is always satisfied and we have zero modes for all 
rational values of $. 
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FIG. 2: (Color online) Energy bands as a function of <1> for t = 0.5 and t = 1.0. 



t=0.5 t=1.0 




FIG. 3: (Color online) A closer look of Figl2] in a weak magnetic field region. The energy levels 
([Mil and (135]) [or (136]) and ([37]) ] are also shown by the red lines. 



B. 1 < t < 2, where zero modes disappear and a gap around E = opens for t > 2^/"?. 

C. t = 2, where one zero mode exists for $ = 0. 

D. t > 2, where no zero modes exist. 



A. < t < 1 

We show two examples of the energy bands as a function of $ in Fig 121 For < t < 1, we 
have zero modes for all rational values of $. As shown in the following, the energy bands in 
a weak magnetic field are well described by the continuum approximation. 

In the continuum approximation, we use the Landau gauge A = (0, Bx, 0) for a magnetic 
field B. Then, substitution Px and Py ^ Py + Bx for ([9]) with p^ = —id^, Py = —idy 

and B = 7r$ (see Appendix gives the equation in a weak magnetic field as 
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where Q± is given by 



v± 

VI 



(27) 



Here, P+ and are those near kQ and fcp , respectively. Since Q± and Py commute each 
other, we can replace Py with a c-number py. Then putting x — x — Py/ir^, we obtain 



(j){x) 



E 



where 



V± 



Q±= \ ^ I , p± = -itp^ ± V4 - f^T^^x. 

From (128!) and (!29|) . the following equation is obtained: 



il){x) 
(p{x) 



ip{x) 
(p{x) 



where 







Therefore, we have 



t^pl + (4 - t2)7r2$2^2 _ _ tH^a^ 



around k^, and 



t^pl + (4 - t^)'K^<l>^x^ + W4: - tH<l>a, 



ipi^x) 

(j){x) 



%Ij{x) 



E' 



%l){x) 
(p{x) 



E 



2 I 



(28) 



(29) 



(30) 



(31) 



(32) 



(33) 



around fcg , where is the 2;-component of the Pauli matrix. Since the equations for if) in 
f l52]) and fl53]) essentially coincide with those for harmonic oscillators, the energy level for ip 
around k^ is given by 



En = ±V27r<l>t(4 - t^)^/^v^, (n = 0,1,2,...), 



(34) 
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FIG. 4: Energy levels around and 



and that around fcg is given by 



E„ = ±V27r$t(4 - t2)^/Vn + 1, (ri = 0,1,2, 



In a similar manner, the energy levels for (f) around kQ and are given by 



En = ±V27r$t(4 - t^Y^Wn + 1, (n = 0, 1, 2, . . .) 



and 



E„ = ±V27r<l't(4 - t^Y^^Vn, (n = 0, 1, 2, . . .) 



(35) 



(36) 



(37) 



respectively. We show the energy levels around k^ and k^ in FigJH As illustrated in FigJSl 
the energy bands for < t < 1 come to be well fitted by (l34l) and ( l35l) [or ( l36l) and ( 1371) ] in 
a weak magnetic field. 



B. 1 <t<2 



For 1 < t < 2, a gap around E = opens for t > 2^/*^. This implies that in a weak 
magnetic field (g ^ 1), a gap around E = opens by a tiny distortion of graphene, t = 1 + e 
(0 < e ^ 1). Since we do not have E = states, the expressions and (!37|) need to be 
modified. We show two examples of energy bands as a function of $ for 1 < t < 2 in Fig|5l 

Let us focus on the gap around E = 0. In order to see how it behaves, we plot the natural 
logarithm of the gap AE around = as a function of 1/$ in Fig|6l From this, we find 
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FIG. 5: (Color online) Energy bands as a function of ^ for t = 1.2 and t = 1.5. 



that it behaves as 



~ exp(-a/<l>). 



(38) 



The values of a are obtained from Fig|6]as a ~ 0.22 and 0.094 for t = 1.2 and 1.5, respec- 
tively. 



30, 



The non-analytic behavior fl38|) can be understood as breaking of supersymmetry 
321] in our model. The operator Q± transforms to (p and vice versa, which is seen from 
( l28l) and ( l29l) . By identifying Q± with generators of supersymmetry, H± in ( l30l) can be 
considered as sypersymmetric Hamiltonians, 7i± = (0 is "boson", and ip is "fermion"). 
Due to the supersymmetry, there is no perturbative (or power-law) correction with respect to 
$ for the E = states. However, tunneling effects break the supersymmetry spontaneously 
and the non-perturbative correction (!38ll appears as a gap around E = 0. 

When two Dirac zero modes at $ = are close to each other in the momentum space, 
namely, t ~ 2, the gap around E = can be estimated by the WKB method. For t ~ 2, the 
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FIG. 6: (Color online) The natural logarithm of the gap around = as a function of 1/$ for 
t = 1.2 and t = 1.5. The slope of the fitting line gives —a. 



two Dirac zero modes at $ = are located at ([7]) with 

cos"^ ^ ~ V2-t = G, (39) 

and for t = 2, they merge into a confluent point ffTTl) . For t ~ 2, it is convenient to expand 
kx and ky around the confluent point (ITTl) instead of or k^: 

kx = kl+P^=7r+P^, ky = k*y+Py=Py, {\Px\,\Py\ <^1). (AO) 

Then T'(fc) in ([5]) is given by 

V{k) = ~2tp^+pl-G^. (41) 

In a weak $, we can use the continuum approximation. We use the Landau gauge 
A = (0, Bx, 0) for a magnetic field B. Then, substitution p^ — > p^, and Py Py + Bx for 
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FIG. 7: (Color online) The asymmetric double-well potential (a) given by (j50p with g = 0.1, and 
(b) given by (j5ip with g = 0.1. Energy levels around g = and q = 1/g are also shown. 



( HTjl with = — "^f^x, = — "ic^y and -B = 7r$ (see Appendix E]) gives the following equation, 



Q 



ipix,y) 



E 



(42) 



where Q is given by 
Q 



V 

V* 



(43) 



In a similar manner as Sec JIV Ai we replace Py with a c-number py and put x ^ x — Py/ir^. 
Then we obtain 



Q 



ip{x) 

(f){x) 



E 



Tpi^x) 
0(x) 



(44) 



where Q is given by 



Q 



P 

V* 



(45) 



Identifying Q with a generator of supersymmetry, we have the supersymmetric Hamiltonian 
H = Q^, which satisfies 



n 



By using the following variable q, 



ip{x) 

0(x) 



E 



2 I ^(^) 



0(x) 



(46) 



G_ _ / 1 



(47) 
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61) can be rewritten as 



with 



1 /7r$ ^ 

g = —\i—. t 



1 



0(g) 



£ 



(48) 



26" V G ' Stt^G 
Therefore, the potential terms for iIj{x) and 0(x) are given by 

1 



(49) 



and 



V_{q) = ^q\l-gqf + gq, 



(50) 



(51) 



respectively (see FigJT]). 

In the leading order of g, the potentials ( 150|) and ( 15T1) are well approximated by the 
harmonic oscillator around q = and q = 1/g. Around q = 0, ( l48l) becomes 



1 rf2 1 2 1 



^(g) 

0(g) 



0(g) 



and around g = 1/g, becomes 



1 
2 



V'(g) 



^(g) 
0(g) 



Therefore, the energy levels for ip around q = 1/g and q = are given by 



S^^ = N+, (iV+ = 0,1,2,...), 



(52) 



(53) 



(54) 



and 



£n.=N. + 1, (iV_ = 0,1,2,...), 



respectively, and those for around q = 1/g and g = are given by 



(55) 



£r,^=N+ + l, (iV+ = 0,1,2,...), 



(56) 
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and 



(iV. = 0,1,2,...), 



(57) 



respectively. 

Now take into account tunneling effects between the energy levels around g = and 
q = 1/ g. The tunneling effects can be estimated by the WKB method presented in Appendix 



D of Ref.j32j. Here we consider only the equation for il){q) since the equation for </)(g) gives 
the same result. The solution of (1521) which vanishes for g — — oo is given by 

(58) 



^{q) = AD, {~V2q 



where v = £ — 1, A is a, constant, and D^, the parabolic cylinder function 331]. The solution 
of (!53|) which vanishes for g — > oo is 



V^(g) = iD,+i \V2{q-l/g) 



(59) 



where A is a constant. We connect these solutions (158!) and (!59|) with that in the forbidden 
region. In the forbidden region, the usual semi-classical expression for the wave function is 
available: 



'4}{q) 



exp (^^ J k{x)da^ + 



^ exp [ f k{x)dx 



(60) 



where k{q) = ^y2{V^{q) — £) with V+(g) in (!50!) . {i = 1,2) are the turning points, 
V+{qi) = £, and Ci {i = 1,2) are constants. Connecting fl58|) with fl60l) . and fl59|) with fl60l) . 
we obtain 



7 



2E-1 



r(i-^)r(-^) = i, 7 



-l/6g2 



1/2 



For £ near zero energy (£^ ^ 1), we have 



2S-1 



r (1-^)^1, r(-£:)^-i. 



thus the solution of ( l6Ti) for <C 1 is obtained as 



1 



7 



87r<l)G' / ' 2 
This implies that the gap around E = is given by 

2G3 



exp 



37r<l> 



(61) 



(62) 



(63) 



(64) 
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FIG. 8: (Color online) a as a function of t, where those obtained from numerical calculations and 
the WKB analysis are shown by points and the line, respectively. 



TABLE I: a for several t obtained from numerical calculations and the WKB analysis. The relative 
discrepancies between them, 5, are also shown. 



t 


a 


awKB 


<5 


1.6 


6.523 X 10^2 


5.368 X 10-2 


0.177 


1.65 


5.198 X 10~2 


4.394 X 10-2 


0.155 


1.7 


4.076 X 10-2 


3.487 X 10-2 


0.145 


1.75 


3.042 X 10-2 


2.653 X 10-2 


0.128 


1.8 


2.140 X 10-2 


1.898 X 10-2 


0.113 


1.85 


1.369 X 10-2 


1.233 X 10-2 


0.0993 


1.9 


7.294 X 10-3 


6.711 X 10-3 


0.0799 


1.95 


2.555 X 10-3 


2.373 X 10-3 


0.0712 



and the exponent a in (l38l) is given by 



a = —G-' 



Stt 



[2 - t)3/2 = awKB, (0 < 2 - t < 1) 



(65) 



Now we compare (!65|) with those obtained numerically for a small $. In FiglHl we show 
a as a function of t, and in Table [H we list them. The relative discrepancy 
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a 



(66) 
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FIG. 9: (Color online) (a)Energy bands as a function of $ for t = 2.0. (b)The three lowest bands 
in £^ > in the log-log scale. 

decreases as t approaches 2. This is because the neglected terms 0{py) in (HTil come to be 
smaller and smaller as t approaches 2. 

C. t = 2 

At t = 2, the two Dirac zero modes at $ = merge into the confluent point (ITT!) . As 
a consequence, a gap around ii^ = in a weak magnetic field makes a transition from an 
exponential (non-perturbative) to a power-law (perturbative) behavior as a function of $. 

In Figini^a), we show the energy bands as a function of $ for t = 2.0. We show the three 
lowest bands in ii^ > in the log-log scale for weak magnetic field in Figj9]^b). We fit our 
data by 

E ~ (67) 

For t = 2.0, we obtain the exponent k as k ~ 0.66, 0.65, and 0.65 for the lowest, the second 
lowest, and the third lowest bands in E > 0, respectively. Thus for t = 2.0 we have a 
power-law behavior E ~ ±$^/^. 

The behavior E ~ ±$2/3 is derived analytically from a particular dispersion relation at 
t = 2 for $ = [25]. For $ = 0, ([6]) and dH]) with G = give 

E = ±y/^^^^, (68) 

which is linear in one direction and quadratic in the other. The exponent k is obtained from 
S{E) ~ $, where S{E) is the area surrounded by an orbit of energy E in the momentum 
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FIG. 10: (Color online) Energy bands as a function of $ for t = 2.5 and t = 3.0. 
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FIG. 11: (Color online) A closer look of FigllOl in a weak magnetic field region. The expressions 
r3|) for n = 0, 1, 2 are also shown by the red lines. 



space |25|. From (!68l) . we have 



(69) 



thus E ~ ±$2/3. 



D. t>2 

For t > 2, we do not have zero modes but have a gap around E = 0. In FigJTOl we show 
two examples of the energy bands as a function of $ for t > 2. 

Let us study behavior of energy bands in a weak magnetic field by the continuum ap- 
proximation. For $ = 0, we expand kx and ky around {k^, ky) = (vr, 0), 

kx = TC+Px, ky=Py, {\Px\,\Py\ <^1), (70) 

then (E]) gives 

E = ±^E^^+2tpl + 2it-2)pl, E, = t-2. (71) 
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FIG. 12: (Color online) (a)Behavior of a gap IS.E around = as a function of t in a weak 
magnetic field. (The nearest bands to -E = show the same behavior as /\E.) (b)Behavior of the 
other energy bands for ii^ ~ as a function of i in a weak magnetic field. 



Thus for small and \py\ (|px|, by| ^ Eg), (17T1) is written as 



E 



(72) 



which is quadratic in both and Py. In the continuum approximation, we put pr^ Px and 
Py ^ Py + Bx with Px = —idx, Py = —idy and B = 7r$. Then we have 



E=±Er. 



l + 2nVi ln + - 



{t - 2)3/2 



(n 



0,1,2, 



(73) 



where we have neglected higher order corrections of O 



V^(n+l/2)$ 
(t-2)3/2 



in the weak magnetic field limit I $ <^ 



The energy bands 
are well approximated by fl73l) . which 



is seen in FigJTTl We note that, for t ~ 2, the neglected higher order corrections of 
Q ^^27rv|(^|n;/2)*~j ^ neglected. However, they become small for t 2, and 

the energy bands in a weak magnetic field are well fitted by 0731) . This result is consis- 
tent with the fact that the honeycomb lattice becomes equivalent to the square lattice for 

t > 2um- 



We summarize our results of this section in Fig 
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V. ^ = EDGE STATES IN ANISOTROPIC HONEYCOMB LATTICE 

In this section, we examine E = edge states. The condition for the existence of zero 
energy edge states in a magnetic field is analytically derived. On the basis of it, it turns out 
that the anisotropy of the hopping integrals induces abrupt changes of the number of zero 
energy edge states, which depend on the shapes of the edges sensitively. 

In order to see this, we focus on the honeycomb lattices with zigzag and bearded edges 
as shown in FigJTSl For these lattices, we have edges along the y-direction. Let us impose 
the periodic boundary condition along the y-direction: 

'4^n,m+2Ly = 'ipn,m, 4'n,m+2Ly = (pn.m, (74) 

where an integer Ly denotes the circumference of the cylinder. Then one can write 

^v™ = exp(4™)v'„, ^„„„ = exp(4™)^„. (75) 

Let us focus on i? = states. In the same manner as Sec JIII Bl for = and ^ = p/q 
with coprime integers p and q, the amplitudes of wavefunctions separated by a distance q 
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are found to satisfy 

\lp(N+l)q+l\ = r\^Nq+l\, \<pNq+l\ = r\(j) (^N +l)q+l\ , 

(/ = 0,l,2,...,g-l), (76) 

with 

_ 2 

From fl76l) . we have 

\ijNq+i\=r''\^il =r^|0;v,+z|. (78) 

The boundary conditions for the zigzag and the bearded edges are given by 

00 = 0, (j)Lx+i = 0, (79) 

respectively. These boundary conditions give 0„ = for all n, thus we consider only ipn in 
the following. 

Suppose that is large enough: ^ q. Then from (!78|) . if r < 1 (r > 1) we have 
E = states on the zigzag edge (bearded edge). For t > 2^/"?, r < 1 is satisfied for all values 
of ky, but r > 1 is not satisfied for any values of ky. Thus we have E = states localized 
on the zigzag edge for all values of ky, but we do not have E = states localized on the 
bearded edge for any values of ky. For t < 2^/"^, E = edge states exist both on the zigzag 
and the bearded edge. The total width dzigzag(t,q) of the region of ky which gives E = 
states on the zigzag edge is given by 

4igzag(t, q)=2(n-2 cos"^ , (80) 

and that for the bearded edge, c?bearded(^, q), is given by 

4earded(t, q) = 2n - 4igzag(i, ?) = 4 COS"^ — . (81) 

For a fixed value of $, (izigzag ((^bearded) increases (decreases) as t increases in the region 
of < t < 2^/'?. At t = 2^/1, 4wag covers the whole region of ky, and c/bearded vanishes. 
For t > 2^^i, we have E = edge states on the zigzag edges for all values of ky. We show 
examples of the energy spectra of honeycomb lattices with zigzag and bearded edges in 
FigsIHandllll 
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cos q 



(77) 



Effects of the anisotropy of the hopping integrals are evident in a weak magnetic field, 
$ ^ 1 (g ^ 1). For t = 1, (izigzag and (ibearded do not depend on the magnetic field, and are 
given by rf^igzag = ^ and (ibearded = respectively. However, for t < 1, d^igzag decreases 
toward and (ibearded increases toward 2tt as q increases. In contrast, for t > 1, (izigzag 
(c^bearded) increases (decreases) as q increases and reaches 2n (0) at g = In 2/ Int. Note that 
even in a small anisotropy, (izigzag and (ibearded change abruptly in a weak magnetic field 

(g»i)- 

Instead of c^zigzag and (ibearded, "we also consider the integrated charge density for i^' = 



edge states 



In = J \Mky)\ dky, (82) 
where the normalization condition is imposed on ipn{ky), 

J^\Mky)\' = l. (83) 

n=0 

To characterize the numbers of the states localized on the edges, we introduce the following 
quantities, 

-^zigzag ^ ^ Iny -^bearded ^ ^ In- (^^) 

n=0 n=Lx—q+l 

For r < 1, only the zigzag edge has E = modes, and A^zigzag is evaluated as 



-^zigzag 



where the domain of the integration is restricted to those ky with r < 1. Here we have used 
the relation derived from flTHl) and fl83l): 



L(L, + l)/gJ-l 

r^' E + 0(r^^^/^) = 1, (86) 

j=Q n=0 

where 0(r^-^"=/'') corrections can be neglected since r < 1 and ^ q. [\x\ denotes the 
integer part of x.) Eq. fl86p is rewritten as 

9-1 

El^„(y|2 = l-r2 + 0(r^^-/''), (87) 

n=0 
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then substituting this into the first equation in (18^ and using (177|) . we obtain Eq. (l85!l . On 
the other hand, for r > 1, we have E = modes only on the bearded edge, and A^bearded is 
given by 

iVbearded = l dky / ^- \dky, (88) 

where the domain of the integration is restricted to those ky with r > 1. Here we have used 
the relation 

|^n(yP = l-r-^ + 0(r-2^-/^), (89) 

which is derived from (1751) and ([HSD in a similar manner as Eq. (1571) . 

Let us now evaluate A^zigzag and A^beardod from fl85|) and fl88l) . For t > 2^/"^, r < 1 is realized 
for all values of ky as shown above. Therefore, 

iVzigzag = ^ ~ t^) ^ / cos{qky)dky 

= 27r ( 1 - 4- 

A^bearded = 0. (90) 

For t < 2^/"^, either r > 1 or r < 1 is realized by a suitable choice of ky. Thus both A^zigzag 
and A^bearded becomc nonzero as. 



2 \ 4 



A^zigzag (^zigzag (1 j_2(7 ) -*2o / „ '"^^ kydky 



2 ( TT - 2 COS ^ — I f 1 - 4-1 + — sin I cos"^ — 

2 y V ^ V V 2 

A'bearded '^bearded „ / TT ~dk,i 

^ .In ^^q2 ( £a 



cos^ I ^ 



^9 / ^9 ~ 

4 cos"^ - - 2t^ sin cos"^ - ) . (91) 



2 V 2 

We also find that for t = 1, A^zigzag and A^bearded are independent of $, 

A'zigzag = 2^3 - A^bearded ^ ^'^ ~ (^^) 

These formulas also show that a small anisotropy induces sudden changes of the edge states 
in a weak magnetic field: When q increases for a fixed t > 1, A^zigzag increases toward 27i and 
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-^Vbearded reaches zero at q = In 2/ Int. And for a fixed t < 1, A^zigzag goes to zero and A^bearded 
increases toward 27r as q increases. Thus in a weak magnetic field (g ^ 1), there are abrupt 
changes in Azigzag and A'bearded at t = 1. 

In FigsJTBl and [T71 we show scaled plots of /„ as a function of n/q. Here we have taken 
Lx = 6g. Because of the normalization condition fl83p . J„ decreases as ~ l/L^ when 
increases. To remove the artificial dependence on L^, we plot g/„ instead_of J„. For t = 1, 

However, for 



g/n's for different $(= p/qYs with the same p fall on a common curve^ 
t > 1, qln comes to be localized on q sites from the zigzag edges as $ (with the same p) 
decreases. In contrast, for t < 1, g/„ comes to be localized on q sites from the bearded 
edges as $ (with the same p) decreases. In Table [TTl we also compare numerical data and 
analytical results presented above. With relative discrepancies less than 10%, they exhibit 
good agreements. 

So far, we have assumed that S> q. When ^ q, the above arguments cannot 
be justified. However, numerical calculations suggest that if L^. ^ 1/$ {pL^ ^ q), the 
particular edge states presented above appear again. In FigUHl we show the energy bands 
for Lx = 50 and $ = 1001/5000(~ 1/5), which are found to be indistinguishable from those 
for $ = 1/5 (FigJHj). It is also found that if the magnetic field decreases and L^. < 1/$ is 
realized, then the energy bands approach those without a magnetic field, which is illustrated 
in FigiH 

VI. SUMMARY AND DISCUSSION 

In this paper, we examined behavior of zero modes and a gap around zero energy for 
a tight-binding model on the anisotropic honeycomb lattice in a magnetic field, whose 
anisotropy is controlled by the hopping parameter t. It was found that zero modes ex- 
ist for all (rational) $ for < t < 1, and a gap around zero energy opens by a tiny 
anisotropy for the graphene in a magnetic field. This is contrasted with the case for the 
square lattice, where zero modes always exist for all (rational) $ when we change the ratio 



of the hopping parameters tx/ty 3J]. For 1 < t < 2, a gap around zero energy in a weak 
magnetic field behaves as a non-perturbative and exponential form as a function of the mag- 
netic field. This non-analytic behavior is naturally explained by tunneling effects between 
energy levels around two Dirac zero modes in the absence of a magnetic field. At t = 2, the 
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FIG. 14: (Color online) Energy spectra of honeycomb lattices with zigzag and bearded edges for 
$ = 1/5 and L^,. = 50. (a) t = 0.9, (b) t = 1.0, (c) t = 2V5(~ 1.15) (d) t = 1.5. The S = edge 
states localized on the zigzag edges are on the blue lines, and those localized on the bearded edges 
are on the red lines. 
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TABLE II: Numerical data, A'^zigzag, -/Vbearded, and analytical results, A^i,,,ardcd- (^) * = 1-0 

(b) t = 1.07 (c) t = 0.97. The relative discrepancies between numerical data and analytical results, 

(^zigzag = |(A^iitLg~ ^zigzag) /-^zigzag I and (5bearded = K^b^ided ~ ^bearded ) /-^bearded I , are also shown. 
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FIG. 15: (Color online) A closer look of FigHlat Ek,{). 

gap around zero energy in a weak magnetic field makes a transition from an exponential 
(non-perturbative) to a power-law (perturbative) behavior as a function of the magnetic 
field. In particular, an explicit form of the gap around zero energy near the transition point 
is obtained by the WKB method. For t > 2, energy bands in a weak magnetic field show 
linear dependence on a magnetic field. 

We also examined edge states with zero energy. The condition for the existence of zero 
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FIG. 16: (Color online) qln as a function of n/q for $ = 1/11, 1/17, and 1/21. (a) t = 1.0, (b) 
t = 1.07, (c) t = 0.97. 
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FIG. 17: (Color online) qln as a function of n/g for $ = 2/11, 2/17, and 2/21. (a) t = 1.0, (b) 
t = 1.07, (c) t = 0.97. 
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FIG. 18: Energy spectra of honeycomb lattices with zigzag and bearded edges for "3> = 1001/5000(~ 
1/5) and = 50. (a) t = 0.9, (b) t = 1.0, (c) t = 2^/^{r^ 1.15) (d) t = 1.5. 

energy edge states in a magnetic field is analytically derived. On the basis of the condition, 
it is found that the anisotropy of the hopping integrals induces abrupt changes of the number 
of zero energy edge states, which depend on the shapes of the edges sensitively. 
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FIG. 19: Energy spectra of honeycomb lattices with zigzag and bearded edges for t = 1.5 and 
= 50. (a) ^> = 1001/5000, (b) ^> = 101/5000, (c) $ = 11/5000, (d) ^> = 1/5000. We also show 
energy bands in the absence of a magnetic field in (e). 



Finally, we would like to discuss possible experimental realization of anisotropy of the 
hopping integrals. Recently, it was experimentally found that a reversible and controlled 



uniaxial strain can be produced in graphene 
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361]. Therefore, we can expect that a small 



anisotropy of the hopping integrals is realized by the uniaxial strain. On the other hand, in 
order to realize a large anisotropy such as t ~ 2, cold atoms in optical honeycomb lattices 
created by laser beams would serve as alternatives 25|]. In these lattices, a large anisotropy 



could be induced and controlled by changing the intensities of the laser fields [37|. In addition, 
by using Raman processes induced by laser fields, effective magnetic fields can be generated 



in the optical lattice 
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{x,y-2b) y 



' -X 

FIG. 20: A honeycomb lattice with a variable lattice spacing. If a = ^/3b, we have the hexagonal 
symmetry for each sublattice. 



Note added. After submission of this paper, we became aware of recent independent work 
which has some overlap with ours for t ~ |4^. We are grateful to G. Montambaux for 
pointing out these papers. 
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APPENDIX A: TIGHT-BINDING MODEL ON A HONEYCOMB LATTICE 
WITH A VARIABLE LATTICE SPACING 

In this appendix, we consider a honeycomb lattice with a variable lattice spacing shown 
in Fig|20l The tight-binding equation in a magnetic field is given by 

Eipix, y) = ty[(j){x + a,y-b) + e^*^«'=0(x + a,y + b)]+ t^(j){x, y), 

E<j){x,y) = tj,[V;(x-a,y + 6)+e-2*'^t(-«)^(a;-a,|/-6)]+t,V'(a;,2/), (Al) 
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where a magnetic flux through a unit hexagon with an area of 5* = 2ab is given by 27r$. Let 
us write 

X = na, y = mb, (A2) 

then (lAip gives 

Eip{na,mb) = ty[0((n + l)a, (m — 1)6) + e^*''*"0((n + l)a, (m + 1)6)] + t^0(na, m6), 

E0(na, mb) = ty[tp{{n - l)a, (m + 1)6) + e-^^'^*^"-^^!!^ - l)a, - 1)^)] + t^H^a, m6).(A3) 

Then, by writing tpn,m = ip{na,mb), (j)n,m = (l){na,mb), and putting t = t^/ty, ty = 1, we 
obtain the tight-binding equation ([T]) from (]A3p . 

When we have lattice spacings a and b in the x and y directions, respectively, as shown 
in Figj20l the momenta and qy in the x and y directions are related to and ky deflned 
in ([3]) as 

k^ = aq^, ky = bqy. (A4) 

Especially, for a = y/Sb, where each sublattice has the hexagonal symmetry, we have 

kx = o-qx, ky = —=qy. (A5) 
V3 

Here, we mention the relation between the magnetic fleld B and the magnetic flux $. 
They are related as 

BS = 2abB = 27r$, (A6) 

that is. 

The tight-binding equation ([T]) corresponds to a = 6 = 1, thus from (]A7[) B and $ are 
related as 

B = 7r$. (A8) 
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